;;----------------------------------------------------------------------
;; plot zonal distribution of ion production 
;;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;load "$KAILIBROOT/lib_radon.ncl"

begin

  plotpath     = "~/fig/"
  plotname     = plotpath+"fig_surf" 
  plotform     = "eps"
  ploteps      = plotname+".eps"
  plotpdf      = plotname+".pdf"
  plottype     = "lat_slice"

  StringFontHeightF  = 0.027 
  cnLineLabelFontHeightF = 0.025 
  tmXBLabelFontHeightF = 0.018
  tmYLLabelFontHeightF = 0.018
  tiMainFontHeightF  = 0.035 
  tiXAxisFontHeightF = 0.035
  tiYAxisFontHeightF = 0.035
  lbLabelFontHeightF = 0.023 

  gsnStringFont = "helvetica-bold"
  gsnStringFont = "helvetica"

  np = 12

  str = (/"","","","","","","","","","","",""/) 

  tiYAxisString = str
  tiYAxisString(0) = "Height above surface (m)" 
  tiYAxisString(3) = "Height above surface (m)" 
  tiYAxisString(6) = "Height above surface (m)" 
  tiYAxisString(9) = "Height above surface (m)" 

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/) 

  abcd = abcd + ") " 

  gsnLeftString = str

  stra = "WCRP1995 " 
  strb = "SW1998 scaled" 
  strc = "Merged " 
 
  gsnLeftString(0) = abcd(0) + stra
  gsnLeftString(1) = abcd(1) + strb
  gsnLeftString(2) = abcd(2) + strc

  gsnLeftString(3) = abcd(3) + stra 
  gsnLeftString(4) = abcd(4) + strb
  gsnLeftString(5) = abcd(5) + strc

  gsnLeftString(6) = abcd(6) + stra 
  gsnLeftString(7) = abcd(7) + strb
  gsnLeftString(8) = abcd(8) + strc 

  gsnLeftString(9)  = abcd(9)  + stra 
  gsnLeftString(10) = abcd(10) + strb
  gsnLeftString(11) = abcd(11) + strc 
 
  gsnCenterString = str

  strx = "                IPRR (cm~S~-3~N~ s~S~-1~N~)"
  stra = strx 
  strb = strx 
  strc = strx 

  gsnCenterString(0) = stra
  gsnCenterString(1) = strb
  gsnCenterString(2) = strc 

  gsnCenterString(3) = stra 
  gsnCenterString(4) = strb
  gsnCenterString(5) = strc 

  gsnCenterString(6) = stra 
  gsnCenterString(7) = strb
  gsnCenterString(8) = strc 

  strx = "         DJF/JJA ratio "
  stra = strx 
  strb = strx 
  strc = strx  
  gsnCenterString(9)  = stra 
  gsnCenterString(10) = strb
  gsnCenterString(11) = strc 

  gsnRightString = str
 
  stra = "ANN" 
  strb = "DJF" 
  strc = "JJA"
  strd = ""

  gsnRightString(0) = stra
  gsnRightString(1) = stra
  gsnRightString(2) = stra 

  gsnRightString(3) = strb 
  gsnRightString(4) = strb
  gsnRightString(5) = strb 

  gsnRightString(6) = strc 
  gsnRightString(7) = strc
  gsnRightString(8) = strc 

  gsnRightString(9)  = strd 
  gsnRightString(10) = strd 
  gsnRightString(11) = strd 

;;----------------------------------------------------------------------
;; load files 
;;----------------------------------------------------------------------

  overland = True
  overland = False
 
  path = "./input/"
  patb = "./input/" 

  run1 = "RA008" 
  run2 = "RA010" 
  run3 = "RA007" 

  filenama = patb+"out_"+run1+"_ANN.nc" 
  filenamb = patb+"out_"+run1+"_DJF.nc" 
  filenamc = patb+"out_"+run1+"_JJA.nc" 

  filenamd = path+"mbase_T63L31.nc" 

  filenamf = patb+"out_"+run2+"_ANN.nc" 
  filenamg = patb+"out_"+run2+"_DJF.nc" 
  filenamh = patb+"out_"+run2+"_JJA.nc" 

  filenami = patb+"out_"+run3+"_ANN.nc" 
  filenamj = patb+"out_"+run3+"_DJF.nc" 
  filenamk = patb+"out_"+run3+"_JJA.nc" 


  fila = addfile(filenama,"r")
  filb = addfile(filenamb,"r")
  filc = addfile(filenamc,"r")
  fild = addfile(filenamd,"r")

  filf = addfile(filenamf,"r")
  filg = addfile(filenamg,"r")
  filh = addfile(filenamh,"r")

  fili = addfile(filenami,"r")
  filj = addfile(filenamj,"r")
  filk = addfile(filenamk,"r")

;;----------------------------------------------------------------------
;; read data 
;;----------------------------------------------------------------------

  rna = fila->RN222    ;; radon  
  rnb = filb->RN222
  rnc = filc->RN222

  rma = fila->RN22290  ;; 90th percentile   
  rmb = filb->RN22290
  rmc = filc->RN22290

  roa = fila->RN22250  ;; 90th percentile   
  rob = filb->RN22250
  roc = filc->RN22250 

  rnf = filf->RN222
  rmf = filf->RN22290 

  rng = filg->RN222
  rmg = filg->RN22290 

  rnh = filh->RN222
  rmh = filh->RN22290 

  rni = fili->RN222
  rmi = fili->RN22290 

  rnj = filj->RN222
  rmj = filj->RN22290 

  rnk = filk->RN222
  rmk = filk->RN22290 

  slm = fild->slm

;;----------------------------------------------------------------------
;; new data structure for zonal average 
;;----------------------------------------------------------------------

  nlev = dimsizes(rna(0,:,0,0)) 

  plev = nlev-1
 
  prna = rna(0,plev,:,:) 
  prnb = rnb(0,plev,:,:) 
  prnc = rnc(0,plev,:,:)

  prnf = rnf(0,plev,:,:) 
  prng = rng(0,plev,:,:) 
  prnh = rnh(0,plev,:,:) 

  prni = rni(0,plev,:,:) 
  prnj = rnj(0,plev,:,:) 
  prnk = rnk(0,plev,:,:) 

  prma = rma(0,plev,:,:)
  prmb = rmb(0,plev,:,:)
  prmc = rmc(0,plev,:,:) 

  prmf = rmf(0,plev,:,:) 
  prmg = rmg(0,plev,:,:) 
  prmh = rmh(0,plev,:,:) 

  prmi = rmi(0,plev,:,:) 
  prmj = rmj(0,plev,:,:) 
  prmk = rmk(0,plev,:,:) 

  proa = roa(0,plev,:,:)
  prob = rob(0,plev,:,:)
  proc = roc(0,plev,:,:) 


  if (overland) then 

  pslm = slm(0,:,:) 

  prna = where(pslm.lt.0.5,-999.,prna)  
  prnb = where(pslm.lt.0.5,-999.,prnb)  
  prnc = where(pslm.lt.0.5,-999.,prnc)  

  prnf = where(pslm.lt.0.5,-999.,prnf)  
  prng = where(pslm.lt.0.5,-999.,prng)  
  prnh = where(pslm.lt.0.5,-999.,prnh)  

  prni = where(pslm.lt.0.5,-999.,prni)  
  prnj = where(pslm.lt.0.5,-999.,prnj)  
  prnk = where(pslm.lt.0.5,-999.,prnk)  

  prma = where(pslm.lt.0.5,-999.,prma)  
  prmb = where(pslm.lt.0.5,-999.,prma)  
  prmc = where(pslm.lt.0.5,-999.,prma)   

  prmf = where(pslm.lt.0.5,-999.,prmf)  
  prmg = where(pslm.lt.0.5,-999.,prmg)  
  prmh = where(pslm.lt.0.5,-999.,prmh)  

  prmi = where(pslm.lt.0.5,-999.,prmi)  
  prmj = where(pslm.lt.0.5,-999.,prmj)  
  prmk = where(pslm.lt.0.5,-999.,prmk)  

  end if 

;;----------------------------------------------------------------------
;; open wks  
;;----------------------------------------------------------------------

  wks = gsn_open_wks(plotform,plotname)               ; open a ps file

;;----------------------------------------------------------------------
;; colormap
;;----------------------------------------------------------------------

  colormap = "amwg"
  colormap = "testcmap"
  colorstart = 50
  colorend   = 190

  colormap = "ViBlGrWhYeOrRe"
  colorstart = 24
  colorend   = 95 

  colormap = "WhViBlGrYeOrReWh"
  colorstart = 5
  colorend   = 80 

  colormap = "WhViBlGrYeOrRe"
  colorstart = 2
  colorend   = 102

 ;gsn_define_colormap(wks,colormap)         ; choose a colormap
  gsn_merge_colormaps(wks,colormap,"gsdtol")         ; choose a colormap

;;----------------------------------------------------------------------
;; plot settings:  
;;----------------------------------------------------------------------

  resp = True
  resp@gsLineThicknessF = 2.0                       ; thickness of lines
  resp@gsLineDashPattern = 2 ;11
  resp@gsLineDashSegLenF = 0.07 ;0.08
  resp@tfPolyDrawOrder = "Draw"

  ;; for all panels 
 
  res                  = True    ; contour mods desired
  res@gsnFrame         = False   ; don't draw yet
  res@gsnDraw          = False   ; don't advance frame yet
  res@gsnMaximize      = False   ; True

  ; map 

  res@gsnMajorLonSpacing = 60.
  res@gsnMajorLatSpacing = 30.
  res@mpFillOn           = False
  res@mpGeophysicalLineColor     = 115
  res@mpGeophysicalLineThicknessF = 0.5
 ;res@mpOutlineDrawOrder = "Draw"

  res@tiMainOn         = False   ; no title
  res@cnFillOn         = True    ; turn on color
  res@cnLinesOn        = False   ; turn off contour lines
  res@cnLineThicknessF = 1.     ; thicker lines
  res@cnLevelSelectionMode  = "ManualLevels"  ; set manual cn levels
  res@cnMinLevelValF   =    0.     ; set min contour level
 ;res@cnMaxLevelValF   = 8000.     ; set max contour level
 ;res@cnLevelSpacingF  = 1000.
  res@cnInfoLabelOn    = False    ; no contour labels
 ;res@cnFillMode       = "RasterFill"  ; turn on raster mode
 ;res@cnFillDrawOrder  = "PreDraw"     ; draw contours before continents

  res@gsnSpreadColors = True    ; use full colormap

  res@gsnSpreadColorStart = colorstart
  res@gsnSpreadColorEnd   = colorend
 
  res@tiYAxisString  = ""  

;;----------------------------------------------------------------------
;; strings 
;;----------------------------------------------------------------------
  res@gsnLeftString    = "" 
  res@gsnCenterString  = "" 
  res@gsnRightString   = ""

  res@gsnLeftStringFontHeightF   = StringFontHeightF 
  res@gsnCenterStringFontHeightF = StringFontHeightF 
  res@gsnRightStringFontHeightF  = StringFontHeightF 

  res@gsnStringFont = gsnStringFont 
 
;;----------------------------------------------------------------------
;; lable bar 
;;----------------------------------------------------------------------
 
  res@lbLabelBarOn       = True       ; no individual label bar
  res@lbLabelStride      = 2          ; every other label bar label
  res@lbOrientation      = "Horizontal" ; vertical label bar
  res@lbLabelFontHeightF = lbLabelFontHeightF 
  res@pmLabelBarWidthF   = 0.7
  res@pmLabelBarHeightF  = 0.08
  res@pmLabelBarOrthogonalPosF =  0.15

  res@trYReverse         = False      ; reverse y axis
 ;res@trYMinF            = 0.         ; set minimum Y-axis value
 ;res@trYMaxF            = 8000.      ; set maximum Y-axis value
 ;res@trXMinF            = 1949.      ; set minimum X-axis value
 ;res@trXMaxF            = 2006.      ; set maximum X-axis value

  res@vpWidthF           = 1.0        ; width of contour plots

;;----------------------------------------------------------------------
;; tickmark 
;;----------------------------------------------------------------------
;;res@tmYLMode = "Explicit"	
;;res@tmYLValues = (/100.,200.,400.,600.,800.,1000.,2000.,3000.,4000./)
;;res@tmYLLabels = (/"100","200","400","600","800","1000","2000","3000","4000"/)
;;res@tmYLMinorValues  = ispan(0,4000,200)

  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; axis  
;;----------------------------------------------------------------------
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF

;;----------------------------------------------------------------------
;; plot array 
;;----------------------------------------------------------------------

  plot = new(np,graphic)
  dum  = new(6 ,graphic)

;;----------------------------------------------------------------------
;; common setting 
;;----------------------------------------------------------------------

  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnMinLevelValF  =   1.0    ; set min contour level
  res@cnMaxLevelValF  =   9.0    ; set max contour level
  res@cnLevelSpacingF =   0.5
  res@lbLabelStride   =   2      ; every other label bar label  


;;----------------------------------------------------------------------
;; ANN panel 1 
;;----------------------------------------------------------------------

  ip = 0 

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prna,res) 

;;----------------------------------------------------------------------
;; ANN panel 2 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnf,res) 

;;----------------------------------------------------------------------
;; ANN panel 3 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prni,res) 

;;----------------------------------------------------------------------
;; DJF panel 1 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnb,res) 

  resp@gsLineColor  = 1                     ; color of lines
  arealons =  40.
  arealone = 180.
  arealats = 50.
  arealate = 80.
  xxpoly = (/arealons,arealone,arealone,arealons,arealons/)
  yypoly = (/arealats,arealats,arealate,arealate,arealats/)
  dum (0) = gsn_add_polyline(wks,plot(ip),xxpoly,yypoly,resp)

;;----------------------------------------------------------------------
;; DJF panel 2 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prng,res) 
  dum (1) = gsn_add_polyline(wks,plot(ip),xxpoly,yypoly,resp)

;;----------------------------------------------------------------------
;; DJF panel 3
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnj,res) 

  resp@gsLineColor  = 1                     ; color of lines

  ; Russia 
  xxpoly = (/arealons,arealone,arealone,arealons,arealons/)
  yypoly = (/arealats,arealats,arealate,arealate,arealats/)
  dum (2) = gsn_add_polyline(wks,plot(ip),xxpoly,yypoly,resp)

  ; Europe
  delete(xxpoly)
  delete(yypoly)
  xxpoly = (/40.,-10.,-10.,40.,40/)
  yypoly = (/75.,75.,40.,40.,75./)
  dum (3) = gsn_add_polyline(wks,plot(ip),xxpoly,yypoly,resp)

  ; China
  delete(xxpoly)
  delete(yypoly)
  xxpoly = (/75,75,120,120/)
  yypoly = (/50,20,20,50/)
  dum (4) = gsn_add_polyline(wks,plot(ip),xxpoly,yypoly,resp)

  ; USA 
  delete(xxpoly)
  delete(yypoly)
  xxpoly = (/240,240,290,290,240/)
  yypoly = (/ 50, 30, 30, 50, 50/)
  dum (5) = gsn_add_polyline(wks,plot(ip),xxpoly,yypoly,resp)
;;----------------------------------------------------------------------
;; JJA panel 1 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnc,res) 

;;----------------------------------------------------------------------
;; JJA panel 2
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnh,res) 

;;----------------------------------------------------------------------
;; JJA panel 3
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnk,res) 


;;----------------------------------------------------------------------
;; DJF/JJA panel 1 
;;----------------------------------------------------------------------


  res@cnMinLevelValF  =   0.     ; set min contour level
  res@cnMaxLevelValF  =   2.     ; set max contour level
  res@cnLevelSpacingF =   0.2

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  pta = prnb
  pta = prnb/prnc

  plot(ip) = gsn_csm_contour_map_ce(wks,pta,res) 

;;----------------------------------------------------------------------
;; DJF/JJA panel 2
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  ptb = prng
  ptb = prng/prnh

  plot(ip) = gsn_csm_contour_map_ce(wks,ptb,res) 

;;----------------------------------------------------------------------
;; DJF/JJA panel 3
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  ptc = prnj
  ptc = prnj/prnk

  plot(ip) = gsn_csm_contour_map_ce(wks,ptc,res) 


;;----------------------------------------------------------------------
;; panel plot 
;;----------------------------------------------------------------------

  resp                    = True   ; panel mods desired
  resp@gsnFrame           = False  ; don't advance frame yet
 ;resp@gsnPanelBottom     = 0.0    ; space for label bar
 ;resp@gsnPanelTop        = 0.5    ; only panel on lower half of page
 ;resp@gsnPanelLabelBar   = True   ; common label bar
 ;resp@pmLabelBarWidthF   = 0.8    ; label bar width
 
  resp@gsnPanelYWhiteSpacePercent = 5
 ;resp@gsnPanelXWhiteSpacePercent = 5
 
 ;gsn_panel(wks,plot(0:np-1),(/4,3/),resp)
  gsn_panel(wks,plot(0:np-1),(/3,3/),resp)

  frame(wks) 


  system("ps2pdf "+ploteps+" "+plotpdf)
  system("mpack -s "+ploteps+" -a "+plotpdf+" k.zhang.iap@gmail.com")

  ;; output plot information 


end
 
          


